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ABSTRACT 

o " 

^vq ■ Aims. We study the hydrodynamic evolution of a non- spherical core-collapse supernova in two spatial dimensions. We begin our 

study from the moment of shock revival - taking into account neutrino heating and cooling, nucleosynthesis, convection, and the 
^ \ standing accretion shock (SASI) instability of the supernova blast - and continue for the first week after the explosion when the 

expanding flow becomes homologous and the ejecta enter the early supernova remnant (SNR) phase. We observe the growth and 
interaction of Richtmyer-Meshkov, Rayleigh-Taylor, and Kelvin-Helmholtz instabilities resulting in an extensive mixing of the heavy 
elements throughout the ejecta. We obtain a series of models at progressively higher resolution and provide a discussion of numerical 
convergence. 

Methods. Diff'erent from previous studies, our computations are performed in a single domain. Periodic mesh mapping is avoided. This 
is made possible by employing cylindrical coordinates, and an adaptive mesh refinement (AMR) strategy in which the computational 
workload (defined as the product of the total number of computational cells and the length of the time step) is monitored and, if 
r^ I necessary, reduced. 

O . Results. Our results are in overall good agreement with the AMR simulations we have reported in the past. We show, however, that 

I . numerical convergence is difficult to achieve, due to the strongly non-linear nature of the problem. Even more importantly, we find 

O ' that our model displays a strong tendency to expand laterally away from the equatorial plane and toward the poles. We demonstrate 

that this expansion is a physical property of the low-mode, SASI instability. Although the SASI operates only within about the first 
C/) , second of the explosion, it leaves behind a large lateral velocity gradient in the post shock layer which affects the evolution for minutes 

C3 . and hours later. This results in a prolate deformation of the ejecta and a fast advection of the highest- velocity ^^Ni-rich material from 

moderate latitudes to the polar regions of our grid within only 300 seconds after core bounce. This effect - if confirmed by 3D 
simulations - might actually be responsible for the global asymmetry of the nickel lines in SN 1987 A. Yet, it also poses difficulties 
for the analysis of 2D SASI-dominated explosions in terms of the maximum nickel velocities, since discretization errors at the poles 
are considered non-negligible. 

Conclusions. The simulations demonstrate that significant radial and lateral motions in the post- shock region, produced by convective 
overturn and the SASI during the early explosion phase, contribute to the evolution for minutes and hours after shock revival. They 
lead to both later clump formation, and a significant prolate deformation of the ejecta which are observed even as late as one week 
after the explosion. This ejecta deformation may be considered final, since the expansion has long become homologous by that time. 
t"-^ ■ As pointed out in the recent analysis by Kjaer et al. , such an ejecta morphology is in good agreement with the observational data 

^^ ' of SN 1987 A. Systematic future studies are needed to investigate how the SASI-induced late-time lateral expansion that we find in 

^^ , this work depends on the dominant mode of the SASI when the early explosion phase ends, and to which extent it is aff'ected by 

the dimensionality of the simulations. The impact on and importance of the SASI for the distribution of iron group nuclei and the 
morphology of the young SNR argues for future three-dimensional explosion and post-explosion studies on singularity-free grids 
that cover the entire sphere. Given the results of our 2D resolution study, present three-dimensional simulations must be regarded as 
k^ \ underresolved, and their conclusions must be verified by a proper numerical convergence analysis in three dimensions. 

^ , Key words, hydrodynamics - instabilities - shock waves - stars: supernovae 
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More than twenty years after the appearance of SN 1987A, Aschenbach et al. 1995; Hughes et al. 2000; Wan g et al 

the landmark event in both modern observational and theo- .Leonard et al.M2006.: .Wang & Wheeler. ■2008.: .Kiaer et al. 

retical supernova (SN) science, the growing observational SN ^^d references therein). 

database continues to challenge supernova modelers. It is cur- 

rently accepted that core-collapse supernovae and their rem- ^j^^ ^„^^^j interpretation of these properties of stellar ex- 

nants are genera ly characterized by large-scale anisotropies j^^^^^^ -^ ^^^^ ^^ ^^ ^^^ ^^^^^ ^^ hydrodynamic instabilities. 

mixmg and smaller-scale clumpmg of material, together with stimulated mainly by SN 1987A, many research groups have 

performed simulations of the late-time (beyond about 1 s after 

Send offprint requests to: T. Plewa core bounce) hydrodynamic evolution of core-collapse SNe in 
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detailed bibliography of ea rly work along these l ines see the first 
paper of the present series. iKifonidis et aDl20Q3h . 

Historically, such simulations can be classified into two main 
types. Models in which a separate multidimensional neutrino 
radiation-hydrodynamic modeling of the supernova explosion 
mechanism is attempted, from which the subsequent late-time 
hydrodynamic evolution is started, belong to the "radiation- 
hydrodynamic explosion initiation" or, in short, the RHD ap- 
proach. In contrast, models in which the detailed physics of the 
"supernova engine" is neglected, and the explosion is artificially 
triggered by some form of manual energy deposition into a stel- 
lar progenitor modelQ may be assigned to the "non-radiation- 
hydrodynamic explosion initiation" or (in short) the non-RHD 
approach. 

The non-RHD approach is the older of the two, and has been 
in exclusive use - in many variations - until the late nineties. It 
continues to be used even today in works where the complicated 
physics that causes the SN explosion cannot be modeled due to 
computational or other constraints, and where the ease in setting 
up a simulation is of primary importance, as e.g. in p erforming 
first studies of three-dimensional phenomena (JHungerford et al. 
l2003l[2005h . in setti ng up several d iff^erent explosions in a num- 
ber of progenitors (^J oggerst et al.l i2QQ9ii) . or in the testing of 
certain hypotheses, such as the possible formation of jets dur- 
ing the explosion and their impact on late-time SN evolution 
(ICouch et al. 2009) . 

Unfortunately, the non-RHD approach is unsatisfactory from 
a theoretical point of view as it neglects physics which is 
known to be important in the supernova context, as e.g. neu- 
trino heating and cooling, and non-radial hydrodynamic in- 
stabilities like the standing accr etion shock instability (or 
SASI, [B londin et al. 2003; Schec k et al.1 l2004l: iLamind l2007l: 
ISch ecket al. 2008; Foglizzo 2009), and multidimensional con- 
vection (Foglizzo et al. 2006; Scheck et al. 2008), which occur 
within the first second after core bounce. 

This has motivated work on the more advanced models of the 
RHD type, which hold the promise to yield a fully consistent hy- 
drodynamic description of the explosion, once the present com- 
putational difficulties w hich are connected to neutrino-driven 
supernova modeling (cf. iBuras et aDl2006bllah have been over- 
come. Su ch RHD models were introduced (in two spatial dimen- 
sions) by IKifonidis et al.l (2000) and improved a nd extended in 
the first two installments of the present series ( Kifonidis et al.l 
I2003LI2006L henceforth Papers I and II, respectively). 

At present, both the RHD and the non-RHD approach make 
use of parametrizations and assumptions. The level of the em- 
ployed approximations is, however, fundamentally diff'erent. In 
the non-RHD approach the entire hydrodynamic structure of a 
model is prescribed. A severe consequence of such "structural 
parametrization", as we will call it in the following, is the ten- 
dency to bias the subsequent evolution with preconceived no- 
tions and assumptions, which are unavoidably introduced in at- 
tempts to set up simple initial data (see the discussion in Papers I 
and II). In contrast, the approximations used in RHD type mod- 
els mainly center on the level of accuracy at which the (com- 
putationally very expensive) neutrino physics is treated. The hy- 
drodynamics is then calculated consistently in response to the 
resulting neutrino heating. The approximations made in RHD 
type models are therefore far less intrusive than those in models 
of the non-RHD type. 



^ usually in one spatial dimension by the employment of a piston 
or a thermal bomb, and sometimes in two dimensions by the use of 
manufactured aspherical initial shock waves 



In Paper I, for instance, a neutrino-hydrodynamics code 
which was based on a simple neutrino light-bulb scheme was 
used to model the evolution encompassing the first second of the 
supernova, given a boundary condition for the dense, contract- 
ing, and neutrino-radiating neutron star core, which was excised 
from the computational domain for reasons of efficiency. 

Paper II improved upon this treatment, by replacing the neu- 
trino light bulb with the more accurate , gray neutrino trans- 
port scheme developed by JScheck et al.l (l2006h . Furthermore, 
the boundary neutrino fluxes were taken to mimic results of su- 
pernova simulations employing full Boltzmann neutrino trans- 
port, but were slightly varied around those values. This led to a 
more accurate description of the ratio of hydrodynamic to heat- 
ing time scales in the post-shock flow than in the models of 
Paper I, and allowed to properly account for large-scale, non- 
radial, low-mode (bi- and quadrupolar) hydrodynamic SASI in- 
stabilities during the first second of the explosion. 

In both Papers I and II the late-time hydrodynamics was sub- 
sequently computed with an Adaptive Mesh Refinement (AMR) 
code. Followed to about six hours after core bounce, the RHD 
models presented in Paper II have demonstrated great potential 
to give good agreement with observations. One particular sim- 
ulation, model b23a, was successful in accounting for the most 
intriguing features of SN 1987A . Based on a rece nt analysis of 
SN 1987A's ejecta morphology Kjaer et al.' ('2010") in fact argue 
that only SASI-dominated explosions (such as model b23a) are 
fully consistent with the observational data , while jet-induced 
explosions, as favored by lWang et al.l (l2002h . are not. 

Of course, the additional physics in such RHD models does 
not come without cost. Temporal and spatial scales need to be 
resolved that are absent in models of the non-RHD type. The 
computational problems encountered in trying to account for this 
physics, along with the need to evolve the models to times suflft- 
ciently late for making comparisons against observations makes 
such detailed modeling daunting, already in two, but even more 
so in three dimensions. 

Th is is exemplified by the difl&culties that iHammer et al.l 
(2009j have experienced in attempting first three-dimensional 
simulations of this type. Their runs were initialized from 3D neu- 
trino radia tion-hyd rodynamic simulations of the explosion per- 
formed by lSchecd (12007,) . which served as data for late-time 3D 
(purely) hydrodynamic calculations. Because of limited com- 
putational resources, the 3D neutrino radiation-hydrodynamic 
models of Scheck (2007) had only followed four nuclear species, 
and lacked a detailed treatment of nucleosynthesis under non- 
NSE conditions. Even with this approximation, his best resolved 
run could only be evolved to a time of 0.5 s after bounce, just 
about half the evolutionary time deemed necessary in 2D mod- 
eling attempts of the explosion phase. The explosion energy in 
his model was therefore far from being saturated. In addition, 
IScheckl (l2007l) had considered boundary neutrino fluxes that 
would have only given a very low-energy explosion. 

In order to still use this model as init ial data for their sub- 
sequent simulations, iHammer et al.l (12009 ) artificially changed 
its energy distribution to match typical SN explosion energies 
of ~ 1 Bethe. Unfortunately, the bias for the later evolution that 
such manipulation introduces is hard to quantify. Such modifi- 
cation of an RHD model might, in fact, be viewed as a regres- 
sion to structu r al para metrization. Moreover, the resolution that 
IHammer et al.l (l2009l) could aff'ord in three dimensions using a 
single-mesh PPM code was naturally substantially inferior to 
that which has been achieved in the best-resolved two dimen- 
sional studies that made use of AMR techniques. Nevertheless, 
their results are interesting, since they indicate that 2D axisym- 
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metric simulations overestimate the drag coefficient of newly 
formed clumps which are enriched in heavy elements, an effect 
which has been predicted for some time (see e.g. Paper I, and 
[Kane et aI]|2QQQ) . While more accurate and better resolved 3D 
studies must and will be performed in the future, there are still a 
lot of unsettled issues even in two spatial dimensions. 

First, few calculations were carried into the phase of homol- 
ogy, and in the cases where this has been done ( Herant & Benzl 
119911: lHerant&Ben3ll992l: iHerant & Wooslevlll994l) the late- 
time models do not include the early explosion stages in a self- 
consistent way. Second, most of the existing models have rela- 
tively low resolution and none of the more recent works known 
to us has discussed numerical model convergence, or provided 
any estimates of accuracy for such important quantities as the 
maximum nickel velocities in the ejecta. As we will show in this 
paper even highly resolved 2D AMR models still show varia- 
tions in the maximum nickel velocities. Third, in models of the 
RHD type the effects on these velocities due to radioactive heat- 
ing by the decay of nickel and cobalt remain to date unquan- 
tified (in models that were started from simple non-RHD ex- 
plosions, Herant & Benz 1992 have noticed a mild acceleration 
due to radioactive heating). Fourth, parameter studies with dif- 
ferent stellar models, to investigate the sensitivity of the results 
on the density and pressure profiles of the progenitor, have not 
yet been performed for models of the RHD type. And fifth, the 
effects of fall-back have not been studied in most simulations, 
since gravity has either been completely neglected or treated in 
the monopole approximation. 

In the present paper we are going to address the first two 
of these points. Our goal is to attain homologous expansion 
in high-resolution mesh-based simulations that are initiated 
from radiation-hydrodynamic models of the early SN explosion 
phase, and to provide a resolution study in order to quantify nu- 
merical uncertainties. Ultimately, we would like to evolve the 
simulations well into the SN remnant phase, in order to compare 
such explosion models - and hence SN theory - with a larger ob- 
servational data base than it has been possible in the past. In the 
process we face some new challenges that such late-time compu- 
tations pose. In particular we aim to resolve, with limited com- 
putational power, and within a reasonable amount of time, all 
crucial physics participating in the evolution. 

Our model s are b ased on the 15 M© blue supergiant star by 
fWoosley et al.' (1988). Within this progenitor is embedded the 
SASI-dominated explosion model b23a that was presented in 
Paper II. To allow for computations beyond the shock breakout 
from the stellar surface, we have added a stellar wind outside 
the progenitor star. We begin our simulations at the final time of 
the explosion model (~ 1 s after core bounce), and continue to a 
final simulation time of 7 days, when the expansion is expected 
to be completely homologous. 

We have structured this paper as follows. In Sect. [2] we shall 
give an overview of our computational approach, and our initial 
and boundary conditions. The results of our simulations will be 
presented in Sect. [3l starting with a discussion of the dynamics 
of spherically symmetric (1-D) models - which serve as a ref- 
erence - and progressing to a high-resolution two-dimensional 
model. In Sect. U we provide a detailed analysis of our mul- 
tidimensional simulations, commenting on aspects such as the 
global lateral expansion of our SASI-dominated models, the re- 
sulting anisotropy of the final metal distribution in the ejecta, 
the extent of the mixing of heavy elements in mass, the evolu- 
tion towards homology, and the numerical convergence in two 
dimensions. Section [5] finally provides a summary of our find- 
ings. 
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2. Computational model 

2.1. Hydrodynamics 

We are numerically solving Euler's equations of com- 
pressible hydrodynamics in conservative form. We use 
the Euleri an version of the P i ecewis e Parabolic Method 
(PPM) by IColella & Woodward! (11984 as implemented in 
Ardent/fla sh, a parallel, blo ck-structured, AMR hydrody- 
namic code (Fryx ell et al.l l200d) . PPM is a high-order shock- 
capturing Godunov-type scheme (lGodunov| [l959), is formally 
second order accurate in space, and uses Stra ng-type directional 
splitting for time advancement (IStran gill 968b to achieve second 
order accuracy in time. We are keeping track of 19 passively ad- 
vected nuclear species by solving a separate advection equation 
for each specie. A list of the species is given in Table [TJ 

The system of hydrodynamic equations is closed using an 
electron-positron equation of state based on table interpolation 
of the Helmholtz free energy (Timmes & Swesty 2000). The 
Helmholtz equation of state (EoS) takes into account contribu- 
tions from electrons and positrons (degenerate and/or relativis- 
tic), radiation pressure, and ions assuming complete ionization. 
The physical limits of this EoS are 10"^^ < p < 10^^ g cm"^, and 
10"^ < T < 10^^ K, appropriate for stellar interiors and the early 
stages of post-explosion supernova remnant evolution. To allow 
for the evolution at very late times, we have extended the origi- 
nal Helmholtz EoS to the low-density, low-temperature regime, 
by blending it smoothly with an ideal gas equation of state. 



2.2. Simulation setup 

Previous simulations of multidimensional core-collapse super- 
novae have been typically performed in spherical coordinates, 
which severely limits the time step due to the small grid spacing 
towards the center, combined with the extreme flow conditions 
prevailing in that region. Another drawback of using spherical 
geometry is the fact that while one may be resolving the interior 
of the star very well, the resolution is quickly lost (especially in 
the angular direction) at large radii. Thus the resolution is also 
highly nonuniform, making convergence analysis difficult. Also, 
while spherical geometry might be optimal at early times of the 
simulation, it is not desirable at later times when the dynamics 
we are most interested in is taking place far away from the cen- 
ter. Yet another disadvantage of using spherical geometry is the 
presence of a grid singularity at the center of the grid, which in 
most cases is avoided by using a computational domain whose 
inner radial boundary is located at a finite radius. We avoid most 
problems associated with a spherical mesh by using a cylindrical 
grid. Also, a cylindrical mesh can be easily combined with adap- 
tive mesh refinement, avoiding cumbersome remeshing/restart 
procedures, and allowing us to take the simulations farther in 
time than we otherwise would have been able to. 

Our two-dimensional cylindrical grid covers the region 
[0, 1.01 X 10^5] cm in radius and [-1.01 x 10^^ 1.01 x 10^^] cm 
in the vertical direction. We impose a reflecting boundary condi- 
tion 3tR = (symmetry axis) and allow for free flow everywhere 
else. The center of the star, where the neutron star is forming, is 
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approximated by a non-moving gravitating sink hole. The hole 
is treated as a point-mass located at the grid center with an initial 
mass of 1.099 M©. It occupies a semi-circular region extending 
up to 3 computational cells from the center. The velocity of the 
material inside the hole is reset to zero after every step. The den- 
sity inside the sink hole is kept constant (^ 10"^ of the density 
of the material just outside the hole) and the mass of material 
accreted into the hole is added to the central point mass (repre- 
senting the neutron star). The gravity of the central mass and the 
self-gravity of the mass residing on the grid are calculated in the 
spherically symmetric (monopole) approximation. 

At the start of our simulation the grid is populated with three 
distinct components. Excluding the inner core occupied by the 
sink hole, the explosion component extends up to a radius of 
r = 1.7 X 10^ cm. (In what follows, the radial distance from the 
grid origin is denoted as r.) The explosion component is the b23a 
model at a time of ^ = 0.92 s after core bounce, that was obtained 
in Paper II. Model b23a in turn is based on the WPE15 Is (180) 
post -collapse model of the 1 5 M© progenitor of IWoosley et al.l 
([T988D presented bv lBruennI (Il993h . 

In an attempt to minimize a possible contribution and accu- 
mulation of numerical discretization errors near the polar axis 
of the explosion model, two 10-degree wide sections near the 
axis of model b23a were overwritten by mirror-reflections of the 
abutting mesh regions (i.e. in these regions the density, for in- 
stance, was initialized according to p{6) = p(20° - 6)). However, 
pronounced axial flows still developed during the subsequent 
evolution. As we will show below, these axial flows are pro- 
moted by the global lateral expansion of that SASI explosion 
model. We recognize that the flow dynamics near the symmetry 
axis is due to a blend of a real physical efl'ect, the fact that the 
degrees of freedom of the flow are restricted, and the likely pres- 
ence of non-negligible discretization errors in these regions (see 
Sect.H. 

The progenitor model is identical to that which we already 
used in Paper II (S. E. Woosley, private communication). The 
model envelope extends up to a radius of ^ 3.9 x 10^^ cm. The 
remaining portion of the computational domain is filled with a 
stellar wind. The wind is spherically symmetric with a velocity 
of 15 km s"\ a temperature of 1 x 10^ K, and a density corre- 
sponding to a constant mass loss rate of 1 x 10"^ M© yr"^. To 
map the explosion and progenitor component to our c omputa- 
tional domain we employ monotonic cubic interpolation (ISteff'enI 
fT99Qh . 



2.3. Computational strategy 

Since Ardent/FLASH uses AMR we are able to study the hy- 
drodynamic evolution with many levels of refinement. In our 
case the mesh is refined if there are density jumps > 0.2, pres- 
sure jumps > 0.4, or if the abundance of nickel exceeds 10% by 
mass. We initially start the simulations with a large number of 
refinement levels and then remove these levels depending on the 
instantaneous number of zones (to limit memory use), and the 
calculated step size (to limit computational time). The latter two 
refinement criteria allow us to control the total computational 
cost of an individual model. 

In Fig. [T] we show the number of mesh cells (bottom curve) 
as a function of the simulation time. The top curve in this figure 
shows the wall-clock time for every double time- step in the sim- 
ulation. We note 4-5 second variations (^15% amplitude) in 
the wall-clock time due to mesh adaptation, which takes place 
every other double time step. Erratic variations of much larger 




log(t) [s] 

Fig. 1. Computational workload in the two-dimensional 15 km 
resolution model of a supernova explosion as a function of the 
logarithm of the simulation time. The number of zones is shown 
with a thick solid line (bottom curve). The simulation begins 
with an eff'ective resolution of 15 km which is gradually de- 
creased with time; at the final time the eff'ective model resolu- 
tion is 393,216 km. The wall-clock time per double time-step 
is shown with a thin line (upper curve). Rapid variations in the 
wall-clock time are partly due to mesh refinement (^15% am- 
plitude change) and code I/O (large amplitude spikes). The sim- 
ulation has been performed at DOE's NERSC on the Franklin 
Cray XT4 supercomputer using 32 cores. 



(^ 20 - 30 s) amplitude correspond to writing output files (ex- 
treme cases are due to I/O congestion on the computing system). 
The general semi-periodic shape of the curves reflects our re- 
finement strategy. During each cycle occupying several hundred 
time steps, the number of computational cells increases. Once 
the maximum allowed number of cells is reached, the finest level 
of refinement is automatically removed and the computation pro- 
ceeds without interruption. The mesh level removal is done grad- 
ually over several time steps to avoid problems with mapping the 
mesh hierarchy to the process space. Occasionally, the flow con- 
ditions may change such that the time step is reduced below a 
certain threshold (e.g. during shock acceleration at log(0 ~ 3.75, 
see also Fig. [2] below) and significantly many more mesh cells 
are removed to limit the computational workload. The total com- 
putational cost of the 15 km model is ^ 1.02 x 10^ steps and 
^ 172 hours on NERSC's Franklin Cray XT4 with 32 cores. 

3. Model results 

Although the main focus of this work is the eff'ects of two- 
dimensional hydrodynamic instabilities, we start with one- 
dimensional models as these are useful for understanding the 
gross overall dynamics of the system, serve as a reference point 
in the development of two-dimensional models, and - very im- 
portantly in the present context - help optimizing computational 
strategies in adaptive simulations. In two spatial dimensions 
fluid instabilities typically grow at contact discontinuities due 
to the preexisting flow perturbations and the interaction with the 



A. Gawryszczak et al.: Non-spherical core collapse supernovae 



\n 15 




11 
log(r) [cm] 



log(t) [s] 



Fig. 2. Supernova shock dynamics in one- and two-dimensional explosion models, (left panel) Evolution of the shock speed (left 
scale) as a function of radius in the one-dimensional (thin) and two-dimensional (medium line) 15 km model. The run of the mass 
distribution (pP) (right scale) in the progenitor star is shown with a thick solid line. Note the tight correlation between the shock 
speed and the mass distribution, (right panel) Temporal evolution of the supernova shock speed (solid lines) and radius (dotted 
lines) in the one-dimensional (thin lines) and two-dimensional (thick lines) 15 km models. Note the rapid shock acceleration after 
the shock reaches the progenitor's outermost layers around log(0 ~ 3.75. 



passing shock wave. This complicates the flow pattern, and usu- 
ally renders the models too costly for making numerous adjust- 
ments to the simulation strategy. 

We expect that our simulations will be largely consistent 
with those reported in Paper II. This is because several essential 
components of the two models, both physics-related and com- 
putational, are common. However, we anticipate to observe also 
some difl'erences due to the highly non-linear nature of the prob- 
lem, the use of a difl'erent grid geometry, the lack of mesh remap- 
ping, and the overall higher mesh resolution. 

3.1. Spherically symmetric model 

We have computed one-dimensional models at efl'ective resolu- 
tions of 240, 120, 60, 30, and 15 km. Figure [3] shows the den- 
sity distribution as a function of the logarithm of the radius in 
these one-dimensional models at ^ = 7 days. As we increase 
the efl'ective resolution from 240 km (dashed line in Fig. [3]) to 
120 km (dotted line), the density increases by ^ 25% at the first 
density maximum (log(r) ^ 13.35) and decreases by ^ 10% at 
the second maximum (log(r) ^ 13.42). Doubling the resolution 
to 60 km (thin solid line) results in a mild density increase at 
both maxima by < 10%. The density at the second maximum 
oscillates with the resolution with a progressively smaller am- 
plitude, while more monotonic convergence is observed at the 
first maximum. Despite the aforementioned variations, the one- 
dimensional density profile seems fairly well established already 
at a resolution of 60 km. 

We begin our simulation at ^ = 0.92 s after core bounce 
with the supernova shock located just past the Si/0 interface at 
log(r) ^ 9.2 cm (see Fig. |4]). Behind the supernova shock one 
can notice a reverse shock at log(r) ^ 8.6 cm that was formed 
when the main shock entered the oxygen core (cf. Paper I). In 
between the two shocks sits a dense shell of metal-rich mate- 
rial which contains most of the ejecta mass. The dense shell is 
moving outward due in part to its momentum, and the combined 
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Fig. 3. Density profiles in log scale as a function of the logarithm 
of the radius in the densest part of the ejecta in one-dimensional 
models at the final time (t = 1 days). Results are shown for 
models with a resolution of 15 (thick solid), 30 (medium solid), 
60 (thin soHd), 120 (dotted), and 240 km (dashed fine). Note 
the gradual decrease in the density variation as the mesh reso- 
lution increases, which is indicative of numerical convergence. 
The density profiles of models with a resolution of at least 60 km 
appear very closely matched. 



action of pressure gradients in the shocked gas, the gravity of 
the central mass, and the self-gravity of the progenitor. The re- 
verse shock present at r = 0.92 s eventually disappears from 
the grid as it passes through the inner region (sink hole), but 
is partially reflected. In addition to the original reverse shock. 
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Fig. 4. Evolution of the density in log scale before shock break- 
out for the one-dimensional explosion model at an effective res- 
olution of 15 km. The initial positions of the Si/0, (C+0)/He, 
and He/H composition interfaces are marked by dotted vertical 
lines. 



two more reverse shocks form during the passage of the su- 
pernova shock through the envelope. The first of those reverse 
shocks originates when the supernova shock enters the helium 
shell, and the second one after the shock enters the hydrogen 
envelope. Those shocks can be seen in Fig. |4] at ^ = 20 s and 
t = 1500 s at radius log(r) ^ 10.25 cm and log(r) ^ 11.68 cm, 
respectively. They form when the leading shock first accelerates 
and then decelerates, which means the freshly shocked material 
is moving slower than the shocked material at earlier time. This 
leads to piling up of the shocked material, and local density and 
pressure increase. If the conditions are right, and the shocked 
material moves supersonically with respect to the dense region, 
this a coustic perturbation becomes a reverse shock (iHerant et al.l 
Il994|) . Similarly to the original reverse shock, the additional re- 
verse shocks are reflected ofl' the central region. 

A shock wave accelerates (decelerates) when it travels 
through a m edium with the density decreasin g faster (slower) 
than oc r"^ (ISedovl [T959I: iHerant et al.1 Il994l) . In other words 
the shock speed depends on the mass overrun by the shock, 
Vsh ^ (pr^)- Due to the density structure of the progenitor, the 
supernova shock propagates through the star in an unsteady fash- 
ion. It is the unsteady propagation of the supernova shock that is 
responsible for Rayleigh-Taylor instabilitie s that develop alon g 
material interfaces dChevalier & Kleinll978L IHerant et al.ll994 . 
The propagation of the blast through the star is illustrated in 
Fig.[2j While moving through the C+0 core, the shock initially 
accelerates to slightly over 20, 000 km s~^ Then it slows down 
in its motion through the helium shell to ^ 10, 000 km s~^ . Upon 
entering the hydrogen envelope, it briefly speeds up, and subse- 
quently decelerates again, to ^ 5500 km s~^ Around log(0 ~ 
3.2 (t ^ 1, 600 s) the shock begins to accelerate in the outer stel- 
lar layers, ultimately reaching velocities of ^ 100, 000 km s"^ 
Once the blast has left the envelope, it gradually slows down 
inside the stellar wind. This last acceleration/deceleration se- 
quence gives rise to the formation of a final reverse shock, and 



is responsible for the growth of a Rayleigh-Taylor instability at 
the interface between the shocked ejecta and the shocked wind 
(see below). The shock is traveling at ^ 32, 500 km s~^ when it 
leaves the computational domain at ^ ^ 3 days. 

3.2. Two-dimensional model 

Although one-dimensional models provide invaluable informa- 
tion about several major characteristics of the exploding star, 
qualitatively new phenomena emerge in multidimensions. These 
are chiefly related to fluid flow instabilities developing at the 
material interfaces that separate various nuclear species. In what 
follows, our intention is to highlight the richness and complex- 
ity of the mixing process rather than to provide a quantitative 
comparison of contributions of the individual participating hy- 
drodynamic instabilities. We feel the value of such a quantitative 
comparison would be extremely limited and could possibly pro- 
duce confusion given the long list of uncertainties involved in 
the physical model (progenitor) or numerics (mesh resolution, 
assumed symmetry). We therefore think it is justified to off'er 
solely a qualitative discussion. 

The main process leading to mixing of diff'erent ele- 
men ts in super novae is the Rayle i gh-Taylo r insta bility (RTI; 
Cha ndrasekha^ 119611: ISharpI Il984l: lYoungsl 119841: ISadotetal.1 
2005). The RTI requires the presence of a material interface sep- 
arating fluids at two difl'erent densities and a sustained acceler- 
ation pointing across the interface from lighter to denser mate- 
rial. This situation naturally arises in thermonuclear supernovae, 
where hot and light ashes of nuclear fuel buoyantly expand into 
unburned mater ial (iNomoto et al.lll976l: iMtiller & Arnet3ll986l: 
lKhokhlovlll995h . In core-collapse supernovae, however, the ac- 
celeration at the material interfaces is due to a positive pres- 
sure gradient that results from the acceleration and decelera- 
tion phases of the supernova shock's motion. As we discussed in 
the previous section, the supernova shock experiences two ma- 
jor deceleration episodes during its propagation inside the enve- 
lope, and one more deceleration phase after it leaves the star and 
sweeps through the wind (see right panel in Fig. [2]). Each decel- 
eration phase creates a positive pressure gradient in layers with 
a negative density gradient, i.e. conditions suitable for RTI. 

Besides the RTI, two more fluid instabilities are at work in 
our case. The Richtmyer-Meshkov instab ility (RMI; iRichtmyeJ 
119601: lMeshkovllT969l: lBrouinettee"2002) occurs at material in- 
terfaces subjected to an impulsive acceleration that (in the core- 
collapse supernova setting) is provided by the supernova shock 
and also a series of reverse shocks sweeping through the ejecta. 
Although the overall radial expansion clearly dominates the ex- 
plosion dynamics, initial flow nonuniformities and later contri- 
butions from the RTI and RMI produce velocity sh ear, naturally 
inducing the Kelvin - Helmholtz instability (KHI: iLambl 11932 : 



ChandrasekhaJl96ll:[Sohn & Hwangll2005l: I Youngs & WilHamsl 



20081) . Both the RTI and RMI, as weU as the KHI, grow from 
small-scale features in the current problem and thus their role 
can be quantified only in well-resolved models. To our knowl- 
edge, no such quantitative comparison of contributions of vari- 
ous instabilities has ever been done. It is also beyond the scope 
of the present study. 

In our two-dimensional models we can only observe the 
combined action of all three instabilities. An assessment of the 
relative contributions of the RMI and the RTI based solely on 
the flow morphology is an extremely daunting task. At the end 
of the linear phase, both instabilities p roduce very similar s truc- 
tures difl'ering only in details (see e.g. lAbarzhi et al.|[2003l) . We 
may expect, however, that difl'erent instabilities will dominate 
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the evolution at various times. The RMI will contribute only af- 
ter a sufficiently aspherical shock interacts with a (radial) den- 
sity gradient. The RTI contribution will be non- vanishing only 
in zones with pressure and density gradients of opposite signs. 
If such conditions can be sustained for a sufficiently long time, 
the RTI with its asymptotic f- growth will (locally) outpace the 
RMI contribution (which scales as ^ f-^). Yet, the present sim- 
ulations - and the RTI growth rates presented in Paper I - in- 
dicate that this is not universally the case. In particular at the 
He/H interface of our progenitor model, RMI-growth dominates. 
This can be understood by considering that the initial growth 
time scales are, trti ~ (length scale/effective acceleration) 2, 
and trmi ~ (length scale)/( velocity difference). The velocity dif- 
ference (postshock vs unshocked envelope) is by far greater than 
the (effective rate of) deceleration at the H/He interface. Also, 
small scale perturbations grow faster in the RMI than in the 
RTI. This latter factor is further amplified by the RMI which 
initially decreases the amplitude of the preexi sting perturbations 
in case of He/H (heavy/light) interfaces (see iBrouilletted 120021 
Sect. 3.1). Both factors result in the RMI dominating over the 
RTI, at least initially. It can also be expected that due to the rel- 
atively weak shear along the interfaces, the KHI growth will be 
modest, will occur at late times, and have less impact than either 
the RTI or RMI on the evolution. 

In our subsequent discussion of the evolution in multidimen- 
sions we will focus on the best resolved model obtained at a 
resolution of 15 km. We will defer commenting on numerical 
convergence to Sect. l4.6l below. Figure [5] shows the density dis- 
tribution at selected times during the first 5 minutes after core 
bounce. At ^ = 4 s (top panel in Fig. O, the central region 
around (z, R) = (0, 0) cm is occupied by the sink hole repre- 
senting the nascent neutron star. The material located right out- 
side the hole and in the equatorial plane is falling back onto 
the neutron star. The supernova shock is highly deformed with 
two low-density bubbles occupying most of the core region, a 
clear reminiscence of the low-mode SASI instability in the b23a 
model (see also Fig. 3(c) of Paper II). The shock front shows 
several kinks (triple points), the most pronounced are visible 
at (z,R) = (-2 X 10^5.8 x 10^), (7.5 x 10^2.6 x 10^), and 
(2 X 10^, 6.6 X 10^) cm. As time progresses, those triple points 
move along the shock surface, and are sites of strong vortic- 
ity deposition at the material interfaces, significantly enhancing 
mixing. 

At ^ = 20 s (Fig. [3 middle panel) the supernova shock has 
already traversed through the (C+0)/He interface and Rayleigh- 
Taylor instabilities are now fully developed in the dense shell 
of material bounded by the supernova shock and the resultant 
reverse shock. The kinks in the supernova shock front are clearly 
visible although their positions are now slightly different. We 
notice clear signs of the Kelvin-Helmholtz instability developing 
at the contact surface associated with the leftmost triple point 
near (-1.1 x 10^^, 2.15 x 10^^). Two large streams of matter form 
a Y-shaped feature along the equator and feed the neutron star. 
The edges of both streams show pronounced KHI decorations. 

During the following few hundred seconds, the supernova 
shock becomes progressively less aspherical while the dense 
central region undergoes dramatic morphological changes. By 
t = 300 s (Fig. O bottom panel), the supernova shock has just 
passed the He/H interface and sweeps through the hydrogen en- 
velope. The slowing shock causes the formation of a "helium 
wall" behind it. Two pronounced "scars" are visible at the wall 
at (-1.4 X 10i\ 1.7 X 10^1) and at (6.5 x 10^^ 1.85 x 10^^). These 
flow features are imprints left by the passing triple points. They 




Fig. 5. Density distribution in log scale in the 15 km resolution 
model, (top) ^ = 4 s; (middle) ^ = 20 s; (bottom) t = 300 s. 
Note the presence of two post-SASI bubbles, wrapped in dense 
metal-rich material, and several triple points (kinks) along the 
forward shock front at early times; the dense shell features well- 
developed RTI decorations at ^ = 20 s. 



will play a crucial role in the deep mixing of the outer layers 
into the central regions with large amounts of hydrogen already 
penetrating into the helium shell. 

At this time we also see that the Rayleigh-Taylor instabilities 
have grown substantially. In two regions, highly decorated RTI 
"fingers" have managed to pass through the layers which are in 
the process of forming the reverse shock at the inner surface of 
the helium wall. The third finger-like feature is located closer to 
the equatorial plane and extends half way through the low den- 
sity core region. The Y-shaped accretion feature, clearly notice- 
ably at earlier times, can now be barely recognized as it becomes 
more extended laterally and suffers from flow instabilities. 

The dense material that penetrates into the layers of the form- 
ing helium wall is of particular interest, as it is composed of fast- 
moving, metal-rich clumps. Heavy elements are also dominant 
in the jet-like polar outflows visible near R = and penetrat- 
ing all the way up to the supernova shock. The origin of these 
features is discussed in Sect. 14.11 

The supernova structure between t = 1500 s and t = 1 d is 
shown in Fig. [6l It can be described as a composite of the su- 
pernova shock and its post-shock region, the mixing layer popu- 
lated with a network of dense, metal-rich structures and inward- 
penetrating hydrogen-rich bubbles, and the central, low density 
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Fig. 6. Density distribution in log scale in the 15 km resolution model, (top left) t = 1500 s; (top right) t = 5000 s; (bottom left) 
^ = 10, 000 s; (bottom right) t = I day. The polar outflows are clearly visible at ^ = 5000 s, afl'ecting model morphology for angles 
< 25° and < 20° from the symmetry axis for negative and positive vertical distances, respectively. The supernova shock position at 
that time is shown with a black solid line. Note that the core morphology is qualitatively similar for t > 5000 s. 



ejecta being swept by the reverse shock formed at the base of the 
helium wall. At ^ = 1500 s (top-left panel in Fig.O, the reverse 
shock is fully formed and is about to overrun the middle RTI 
feature. Hydrogen bubbles are now more clearly visible, espe- 
cially the one in the upper left section of the ejecta. Most of the 
fast moving metal-rich material is trapped inside the helium wall 
although some material (z,R) = (3.2 x 10^\ 5.8 x 10^^) appears 
successful in penetrating all the way into the hydrogen envelope. 

By t = 5000 s (top-right panel in Fig.O the supernova shock 
is accelerating in its motion through the outer parts of the stel- 
lar envelope. It breaks through the stellar surface at r ^ 6000 s. 
Substantially more metal-rich clumps are now moving through 
the hydrogen envelope and both hydrogen bubbles are fully de- 
veloped. There are signs of a mild RTI developing at the outer 
surface of the helium wall, best visible close to the equator. The 
reverse shock is sweeping through the low density core overrun- 
ning the network of the RTI-fragmented, metal-rich shell in the 
process. The central part of the ejecta undergoes relatively little 
changes later. Most metal-rich features visible ait = 5000 s can 
be clearly identified also at r = 10, 000 s, and t = 1 day (bottom- 
left and bottom-right panel in Fig.O respectively). The only ma- 
jor change is due to the passage of the reverse shock through the 
central region of the ejecta. This event does not aff'ect the overall 
structure of the core, however. 

Although at these late times the ejecta morphology appears 
well-established, it is the outer regions of the young supernova 
remnant where a significant RTI finally comes into play. This 
is due to the final deceleration of the supernova shock inside 
the stellar wind: the less dense shocked wind material is now 
moving slower than the shocked and denser outer stellar layers. 
The interface separating the shocked stel lar envelope from the 
shocked stellar wind is RTI-unstable (see IChevalier et al]| 19921: 
iNymark et al.ll2006l and references therein) and numerous RTI 
fingers are formed in the process (Fig. [7]). 

Figure [8] shows the angle- averaged radial structure of the 




z(xl0^12cm) 



Fig. 7. Density distribution in log scale in the 15 km resolution 
model at ^ = 10, 000 s including the supernova shock. Note the 
Rayleigh-Taylor instability developing at the interface between 
the shocked stellar ejecta and the shocked wind at a distance 
r ^ 2.4 X 10^^ cm from the explosion center. Note that at this 
intermediate time, the model morphology at large radii is clearly 
aff'ected by the polar outflows for angles < 35° and < 15° from 
the symmetry axis for negative and positive vertical distances, 
respectively. 



young supernova remnant at t = 12 hours. The outermost RTI- 
mixing region is clearly visible in the density, velocity and pres- 
sure profiles for r = [1.3 x 10^"^, 1.7 x 10^"^] cm. Note that at 
this time the metal-rich core occupies only ^ 10% of the inner- 
most region (r < 2 x 10^^ cm). The ejecta are now approaching 
homologous expansion as indicated by the low ratio of the inter- 
nal to the total energy. The only exception is the freshly shocked 
RTI-unstable outermost region. The relatively large internal en- 
ergy content of the core (^ 10% at r = 1 x 10^^ cm, and a few 
percent at the core's edge) is due to the recent passage of the re- 
verse shock; the core continues to expand ballistically and will 
quickly cool adiabatically. 
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Fig. 8. Angle-averaged structure of the 15 km model at ^ = 12 
hours. Shown is the density (thick solid), pressure (medium 
solid), radial velocity (thin solid), hydrogen mass fraction (thin 
dotted), nickel mass fraction (thick dotted), and ratio of internal 
to total energy (thin solid). The RTI-unstable interface between 
the shocked ejecta and the shocked stellar wind is located around 
r ^ 1.5 X 10^^ cm. The metal-rich core extends from the center 
up to r ^ 1.5 X 10^^ cm. Note that additional scaling factors are 
used for the pressure, nickel abundance, and the energy ratio. 
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Fig. 9. Temporal evolution of integrated quantities in the two- 
dimensional 15 km resolution model, (thick solid) explosion en- 
ergy; (medium solid) total energy (sum of the internal and ki- 
netic energy); (thin solid) gravitational energy; (thin dotted) in- 
ternal energy; (medium dotted) kinetic energy; (thick dotted) in- 
tegrated 0-component of the vorticity, cl)^. The data are shown 
up to the moment when the supernova shock leaves the compu- 
tational domain atr^2.6xl0^s. 



The low overall internal energy content is confirmed in Fig. [9] 
in which we show the evolution of the model's energies and the 
integrated vorticity up to the moment when the shock leaves the 
computational domain. The internal and kinetic energy vary to- 
gether in sync, reflecting the variations in the shock speed (cf. 
Fig.O. They are both ^10^^ ergs as long as the shock remains 
inside the star. The internal energy rapidly decreases after the 
shock breakout through the stellar photosphere, and amounts to 
^ 1.5% of the total energy at ^ ^ 23 hours, when the shock leaves 
the grid. It is interesting to note the large variations of the vor- 
ticity during the shock breakout. These variations are due to the 
rapid increase in the shock speed and slight global asphericity 
of the shock front. The vorticity evolution prior to shock break- 
out is consistent with that of the models presented in Paper II. A 
direct comparison of our Fig.[9]with Fig. 5 of Paper II is not pos- 
sible, though, as the latter displays the evolution of the vorticity 
in inner layers of the ejecta (i.e. up to the He/H interface) for 
the early phase after shock breakout. In passing we note that a 
realistic model of shock brea kout must include radiation efl'ects 
(ICalzavara & Matzne3l2004 and references therein), which are 
not considered in this work. We expect, however, that these ef- 
fects will not afl'ect mixing of the metal-rich core nor the struc- 
ture of the young supernova remnant emerging within the first 
few days after explosion. 

The density distribution in the model at the final time, ^ = 7 
days, is shown in Fig. [TOl Except for the apparent density in- 
crease in the central region due to the reverse shock, the over- 
all structure of the metal-rich part of the ejecta is quite simi- 
lar to that obtained a few hours after the core bounce, and the 
major high- velocity features can be easily identified already at 
t = 5000 s (see Fig. O. We also note that the origins of some 
ejecta features, in particular those of the hydrogen-rich bubbles. 




Fig. 10. Density distribution in log scale in the 15 km resolution 
model at ^ = 7 days. 



can be tracked all the way back to the first few minutes after the 
core bounce. 

The most conspicuous feature visible in Fig. [TOl though, is 
the pronounced anisotropy of the ejecta. These have approxi- 
mately the form a prolate ellipsoid with a major to minor axis 
ratio around 1.6. This was already pointed out in Paper II based 
on simulations that were evolved to a time of 20 000 seconds 
after core bounce. Figure [TOl demonstrates that this result holds 
even for the evolution into the homologous phase. In fact, Fig.fTOl 
might be regarded as giving an accurate picture of the ejecta 
morphology even beyond that phase since the ensuing expansion 
is expected to be self- similar, i.e. the morphology will remain 
frozen in the flow. 
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Fig. 11. Angular structure of the 15 km resolution model, (left panel) radial averages of density (thick solid), pressure (thin solid), 
radial velocity (thin dotted), and lateral ^-velocity component (thick dotted) at the initial time, t = 0.92 s; (right panel) radial average 
of density at 0.92 s (thin solid), 1500 s (thin dotted), 5000 s (medium dotted), 10, 000 s (thick dotted), 12 hours (medium solid), 
and 7 days (thick solid). Radial averages in the right panel include only data for a distance < 1.007 x 10^^ cm from the origin. 
Note the gradual evacuation of the equatorial regions of the model and the pile-up of material near the poles (right panel), and the 
strongly anisotropic lateral velocity distribution in the explosion model (left panel, thick dotted line) characteristic of a low-mode 
SASI convective pattern. 



4. Discussion 

4.1. The SASI-induced lateral expansion: Flow dynamics 

A very important effect that we find in the present simulations, 
whose consequences have not been recognized previously (al- 
though they were present in the simulations of Paper II), is the 
presence of a large lateral velocity gradient in our initial SASI 
explosion model (cf. the left panel of Fig. [TT]). Due to this gra- 
dient the model shows a strong lateral expansion, away from the 
equator and toward the poles which leads to an accumulation of 
material in the polar regions, within minutes to hours after core 
bounce (see the right panel of Fig.fTH. 

Figure [121- which depicts the angular distribution of several 
radially-averaged quantities in the young supernova remnant at 
t = 1 days - demonstrates that this anisotropy is preserved until 
the end of our simulations. Note that the density (shown with 
a thick solid line in Fig. [121 and in the left panel of Fig. [TT]) is 
highest inside the polar regions. Note also that the initial angular 
distribution of the density is essentially uniform (shown with a 
thin solid line in the right panel of Fig.[TT1). 

The global lateral expansion contributes decisively to shape 
the ejecta into the final form visible in Fig.[TOl It is also essential 
for the formation of the bulges of material that are visible at the 
poles of the exploding star as early as ^ = 300 s into the evolu- 
tion (bottom panel in Fig.O. These axial flows are the result of a 
complex interaction between the strong SASI lateral expansion, 
and the fact that reflecting boundaries are used at the symme- 
try axis which restrict the degrees of freedom of the flow. Being 
continuously fed by the lateral motion, and encountering the im- 
penetrable boundary at the axis, the path of least resistance for 
the fluid is to expand along the poles. 

In the past, these axial flows have been regarded to be mainly 
the result of discretization errors near the axial coordinate singu- 
larity, and conical sections near the poles were used to exclude 
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Fig. 12. Angular structure of the 15 km resolution model at the fi- 
nal time. Radial averages are shown for the density (thick solid), 
pressure (medium solid), radial velocity (thin solid), and mass 
fraction of hydrogen, oxygen, and nickel (dashed, thin dotted, 
and thick dotted line, respectively). Radial averages include only 
data for distances < 2.4 x 10^^ cm from the origin. 



these regions from further analysis. However, the pivotal role 
of the SASI-induced lateral expansion, indicates that this inter- 
pretation does not account for the true contribution of the vari- 
ous involved eff'ects. Without the strong continuous feeding by 
the global lateral flow, the axial outflows cannot become that 
strong. This is indicated by the fact that the SASI-dominated 
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explosions which are considered in both the present paper and 
Paper II show very pronounced axial flows, while the models 
presented in Paper I, where the SASI was absent, showed much 
weaker axial features. Consistent with this observation is the fact 
that in the limit of a completely vanishing lateral flow, as it is 
for instance the case in 2D simulations of Sedov blast waves on 
spherical grids, such axial outflows are not present. 



4.2. The SASI-induced lateral expansion: Heavy elements 
redistribution 

Another very important efl'ect of the SASI-generated lateral ex- 
pansion is that it causes a lateral redistribution of material en- 
riched in heavy elements. Figure [13] shows the evolution of the 
angular distribution of the nickel mass in our 15 km resolu- 
tion 2D model. The initial SASI explosion model contains two 
nickel-rich sectors centered around ^ ^ 45° and 6 ^ 90° (thin 
dotted line in the left panel of Fig. (TS]). During the subsequent 
evolution, the SASI-generated lateral velocity field of this model 
(cf. Fig. [TT]) causes a steady drift of the nickel toward the polar 
regions. 

Of the two nickel concentrations, the material initially lo- 
cated at ^ ^ 45° is actually the fastest moving nickel-rich ma- 
terial on our grid, with velocities approaching 4000 km s~^ (see 
the right panel of Fig. (TS]). This material is advected over half a 
quadrant of the computational domain and to the north pole in 
less than 300 seconds after core bounce (thin solid line in the left 
panel ofFig.[T3]). 

The evolution of the nickel-rich sector initially located near 
the equator is somewhat more complex as it eventually dis- 
tributes nickel in a broad region between ^ 135° and ^ 
165°. There is a relatively well-defined nickel-rich "clump" (ac- 
tually a torus-like structure, given that our simulations make 
use of axisymmetry) at ^ ^ 145° moving with a velocity of 
;^ 2500 km s-^ 

In addition to becoming enriched with nickel, the polar re- 
gions are relatively rich in oxygen (thin dotted line in Fig. [T2l) 
but devoid of hydrogen (dashed line in Fig. [121). This lateral ad- 
vection of heavy elements poses difficulties for the analysis of 
our 2D simulations, as it transports the nickel clumps over time 
into regions of the computational domain where discretization 
errors are expected to be non-negligible and the accuracy of the 
solution is largely unknown. There is significant ambiguity in 
what should and what shouldn't be accounted for in these re- 
gions when computing the maximum nickel velocities. Given 
that it is not just numerical errors which are contributing, but a 
true physical eff'ect (i.e. the global lateral expansion) is involved 
too, it is for instance not at all clear whether the use of an ex- 
clusion cone around the symmetry axis is a justified and viable 
approach, or how large that exclusion cone should be made. This 
ambiguity is impossible to overcome in axisymmetric models. It 
is of a particular consequence since in our highest resolution 2D 
model it aff'ects the highest velocity nickel clumps present in our 
2D simulations which are of importance for explaining the data 
ofSN1987A. 

To illustrate this further we show in Fig. [141 the mass distri- 
bution of the nickel in velocity space for our 15 km resolution 
2D model at ^ = 7 days, as obtained by employing exclusion 
cones of diff'erent widths. The maximum nickel velocity as de- 
termined without using any exclusion cone at all (A^ = 0°), is 
close to 4000 km s~^ This value decreases gradually to 3350 
km s~\ 3250 km s"^ and 3200 km s"^ if cones with a width A^ 
of 10°, 20°, and 30° are used, respectively. The error in deter- 



mining this velocity could thus be as high as 25%, depending on 
what one considers to be a true feature of the solution. 
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Fig. 14. Mass distribution of the nickel in radial velocity space 
in the two-dimensional 15 km resolution model at ^ = 7 days. 
Data for diff'erent widths of the polar exclusion cone A^ (solid: 
0°, 10°, 20°, and 30°; dashed: 22.5°) are shown. 



4.3. Extent of the mixing of different elements 

In evaluating the extent of the mixing of diff'erent elements in 
mass we again show results considering both the entire computa- 
tional domain, and a 135° wide wedge centered around the equa- 
tor (z = cm) (thereby excluding a cone with width A^ = 22.5° 
around the symmetry axis). The latter approach is similar to what 
has been done in previous work (though the conical section used, 
e.g., in Paper II had an extent of only 15°). 

Figure \T5\ compares the amount of mixing in the 15 km 
model as inferred from the entire domain (left panel) to the mix- 
ing in the region limited to the equatorial wedge (right panel). As 
we can see, essentially all material mixed beyond Mr ^ 9 Mq is 
associated with the polar regions of the model. Slightly increas- 
ing the extent of the exclusion cone, i.e. making the equatorial 
wedge still narrower, does not change the above picture signifi- 
cantly. 

Comparing the species distribution as a function of mass at 
^ = 10, 000 s (right panel in Fig. \T5^ with Fig. 6 of Paper II, we 
observe that at this intermediate time the extent of mixing in our 
present high-resolution two-dimensional simulation matches the 
amount reported in Paper II very well. For example the average 
position of the He/H interface given in Paper II is ^ 5 M©, es- 
sentially identical to that obtained in our 15 km model. In both 
studies H is mixed down to 1.3 M© while ^^Ni is mixed out to 
^ 9 Mq. The agreement is thus excellent. 

By comparing the species distribution at r = 10, 000 s (right 
panel in Fig. [151) to that at ^ = 7 days (Fig. [161), we can, moreover, 
see that most of the mixing takes place within the first 2 or 3 
hours of the supernova explosion. This observation is consistent 
with the ratio of the internal to total energy decreasing to ^ 0.15 
by t = 10, 000 s (see Fig. [island the discussion in Sect. [13). At 
that time, lateral pressure gradients are very small, and the fluid 
elements move nearly ballistically. 

As it has been discussed extensively in Paper II, both the 
inward mixing of hydrogen and the outward mixing of nickel 
are crucial for an understanding of SN 1987A's data. The in- 
ward mixing of hydrogen down to 1.3 M© is required for ex- 
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Fig. 13. Distribution of nickel in the 15 km resolution model, (left panel) The angular distribution of the nickel mass is shown at 
0.92 s (thin dotted), 20 s (thick dotted), 300 s (thin solid), 5000 s (medium solid), and ^ = 7 days (thick solid). Note the fast 
advection of nickel to the north pole (6 ^ 0°). (right panel) Angular distribution of the nickel mass as a function of distance from 
the origin (right axis) and the radial velocity (left axis) at ^ = 7 days. The nickel mass and average radial velocity are calculated 
using bins with Ar = 2.4 x 10^^ cm and A^ = 3°. The most abundant bins are marked, with a nickel mass of at least 1 x 10""^ M© 
(solid circles) and 1 x 10"^ M© (open circles). The average radial velocity for bins located at the same radial distance is shown with 
dotted lines. Note again the large amounts of high velocity nickel located near the symmetry axis at ^ = 0° and the less massive 
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plaining the broad peak of the light curve and the very low min- 
imum hydrogen velocities of this supernova, which according to 
iKoz ma & Fransson (1998) were < 700 km s~^ FigurefTTl which 
displays mass distributions in velocity space for the different 
nuclear species at the final simulation time, demonstrates that 
this is excellently reproduced by our models. Again an exclu- 
sion cone around the symmetry axis with A^ = 22.5° has been 
used for producing this figure. The outward mixing of nickel 
to a mass coordinate of 9 M© is tantamount to the presence of 
nickel-rich material on the grid with velocities > 3200 km s~^ 
We will return to this point in Sect. 14.61 where we will discuss 



the dependence of the nickel velocity distribution on the numer- 
ical resolution. 



4.4. Evolution towards homology 

The evolution of the system towards homology is character- 
ized by the ratio of internal to total energy. As we discussed in 
Sect. 13.21 the internal energy of the exploding supernova varies 
during the early stages of shock propagation through the enve- 
lope. It begins to decrease steadily for t > 1500 s (see Fig.O, 
when the supernova shock starts its final acceleration phase prior 
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Fig. 17. Mass distributions in radial velocity space in the two- 
dimensional 15 km resolution model at ^ = 7 days. The distribu- 
tion of mass is shown for hydrogen (thin solid), ^He (thin solid 
with open circles), ^^O (dashed), ^^Si (dash-dotted), "^"^Ti (dot- 
ted), and ^^Ni (thick solid). Only data for the equatorial section 
of the model are shown (polar regions are excluded). 



to breaking out through the stellar surface. At this time the in- 
ternal energy amounts to about half of the total energy of the 
supernova and exhibits substantial variations in angle (medium 
solid line in Fig.fTSl). 

These are largely smoothed out hy t = 5000 s (thin dotted 
line in Fig. [18]). The average internal energy content drops to 
^ 15% by t = 10,000 s, and well below 5% by ^ = 12 hours. 
It is unclear whether this remaining energy is sufficient to mod- 
ify the expansion velocities of our 15 km model. Using a polar 
exclusion cone of 22.5° width around the poles, we observe a 
relatively modest increase of the maximum nickel velocity on 
the grid by about 12%, from ^ 2850 km s"^ at t = 10,000 s 
(not shown) to ^ 3200 km s"^ at ^ = 7 days (Fig.fTTl). However, 
and as explained in Sect. 14.21 at these late times the nickel-rich 
material in question is located so close to the north pole, that the 
determination of the maximum nickel velocity is affected by the 
ambiguity in choosing the width of the exclusion cone (see also 
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Fig. 18. Evolution of the internal to total (internal plus kinetic) 
energy ratio in the 15 km resolution model. Radial averages of 
the energy ratio are shown at 0.92 s (thin solid), 1500 s (medium 
solid), 5000 s (thin dotted), 10, 000 s (medium dotted), 12 hours 
(thick solid), and 7 days (thick dotted). Radial averages include 
only data for distances < 2.517 x 10^^ cm from the origin. The 
apparent excess heat content of the explosion model (thin solid 
line) near the equator is due to low fluid velocities in a funnel 
separating two rolls (bubbles) which were created by the SASI 
in the initial phases of the explosion. 



Sect. 14.61) . With the morphology of the flow remaining essen- 
tially unchanged for t > 5000 s, it is, however, safe to state that, 
for the employed blue supergiant progenitor model, the ejecta 
become homologous during the first day of the evolution. 

4.5. Comparison to one-dimensional models 

As a final illustration of the importance of the multidimensional 
eff'ects in our models we contrast angle- averaged radial profiles 
obtained from our 15 km resolution two-dimensional model (left 
panel in Fig. [191) with plots of our corresponding spherically 
symmetric (ID) model with the same resolution and at the final 
time (right panel in Fig.[T9b. 

Both models appear to produce explosions of comparable en- 
ergies as indicated by the radial velocity profiles. Both models 
are also quite similar for radii > 2 x 10^^ cm, i.e. in regions unaf- 
fected by mixing in multidimensions. However, at smaller radii, 
we observe an extensive mixing of nuclear species. For example, 
substantial amounts of oxygen are present at r ^ 1.5 x 10^"^ cm, 
and hydrogen is mixed down to r ^ 3 x 10^^ cm in the multi- 
dimensional model. The density jump associated with the H/He 
interface, which is clearly visible at r ^ 1.15 x 10^^ cm in 1-D, 
is completely smeared in multidimensions. The only source of 
mixing in 1-D is numerical diff'usion, while hydrodynamic in- 
stabilities - which are partly caused by the initial shock nonuni- 
formity - contribute to mixing in multidimensions. 

Note the presence of regions of opposite density and pres- 
sure gradients, as required by the Rayleigh-Taylor instability, 
in the spherically symmetric model atr^6xl0^^cm and 
r ^ 1.1 X 10^^ cm. The former region is associated with the 
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Fig. 19. Structure of the exploding supernova in the 15 km resolution models at the final time of t = 1 days, (left panel) Two- 
dimensional model, (right panel) One-dimensional model. In both panels we show the distributions of density (thick solid), pressure 
(medium solid), radial velocity (thin), hydrogen mass fraction (dashed), ^^O mass fraction (thin dotted), and ^^Ni mass fraction (thick 
dotted line). In addition, the ratio of internal to total energy in the two-dimensional model is shown in the left panel (dash-dotted). 



outer edge of the nickel-rich region, while the latter is connected 
to the H/He interface; both regions experienced vigorous mixing 
in multidimensions, as expected. 



4.6. Numerical convergence in multidimensions 

Our final concern is mesh convergence (see Sect. 13.11 for a dis- 
cussion of the one-dimensional case). The quantity which is of 
largest importance for comparison to observations and which 
shows the most sensitive dependence on the mesh resolution of 
our 2D models is the nickel distribution. In contrast, the distribu- 
tion of e.g. hydrogen is very similar in all of our 2D models. We 
will therefore focus on the nickel distribution in what follows. 

It is important to note that for the present problem proper 
convergence analysis is made very difficult by several eff'ects. 
First, no analytic solution for the problem is known, making it 
impossible to determine the actual solution error. Second, the 
problem is highly non-linear. Thus, even if one would obtain an 
estimate of the instantaneous solution error by Richardson ex- 
trapolation, the contribution and interaction of higher order error 
terms, which are not captured by this technique, must be con- 
sidered highly likely. An estimate of the instantaneous solution 
error is, moreover, not relevant as we are primarily interested in 
the accuracy of the solution at late times, when the error has ac- 
cumulated. And third, in the context of the present axisymmetric 
simulations, significant difficulties and ambiguities for the anal- 
ysis of the simulations are encountered near the poles. 

In Fig. [141 we have already demonstrated that the high- 
velocity tail of the nickel distribution of our 15 km resolution 
model at ^ = 7 days depends sensitively on the width of the ex- 
clusion cone that is used around the symmetry axis. Figure [2Ql 
shows that the same is true in our lower-resolved 30 km and 
60 km models. Particular noteworthy is that the maximum nickel 
velocity as obtained from these figures depends on both the 
width of the exclusion cone, AG, and on the numerical resolu- 
tion of the simulation. 

Table [2] summarizes this dependence for convenience. It 
demonstrates that for no value of AG, convergence in the maxi- 
mum nickel velocities is observed with increasing resolution. In 



Table 2. Dependence of the maximum nickel velocities (in 
kms~^) in our two-dimensional models with diff'erent spatial 
resolution on the width of the polar exclusion cone A^. 



Resolution 



0° 



AG 
10° 20° 



30° 



15 km 3950 
30 km 3800 
60 km 3600 



3350 
3800 
3425 



3250 3200 
3450 3400 
3350 3350 



fact, the obtained maximum nickel velocities are not even mono- 
tonic if an exclusion cone with a finite width is used. Only for 
the case where no exclusion cone at all is employed (A^ = 0°) 
do we observe monotonically increasing maximum nickel veloc- 
ities with increasing resolution, which, however, do not show a 
clear sign of convergence. 

From Table [2l we can only infer a lower limit for the max- 
imum nickel velocities at ^ = 7 days of > 3200 km s~^ It is 
also noteworthy that for a 240 km resolution model, GuzmanI 
(2009) reported a qualitatively diff'erent nickel distribution, with 
no nickel moving faster than ^2100 km s~^ In his study, signif- 
icant amounts of nickel moving with velocities > 3000 km s~^ 
are observed for the first time at a resolution of 120 km, which 
is consistent with what we find here. The maximum velocity of 
Fe-group elements reported in Paper II, where an exclusion cone 
of A^ = 15° was used, is ^ 3300 km s"^(see Fig. 8 of that work). 
This is in good agreement with the above results, although mesh 
remapping, and a diff'erent mesh geometry were used in this ear- 
lier work, and the evolution was followed to only 20,000 seconds 
after core bounce. 

From Figs. [14] and [20l it is furthermore obvious that also 
the amount of nickel in the high- velocity tail shows no signs 
of convergence. Even between the better resolved models it can 
diff'er by up to a factor of three, when no exclusion cone is 
used to analyze the simulations. Using exclusion cones of fi- 
nite width, the diff'erences become even larger. In contrast, the 
mass of heavy elements moving at 1000 - 2400 km s"^ is quite 
similar in the models presented here: the mass of nickel expand- 
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Fig. 20. Same as Fig.[T4lbut for our two-dimensional models with 30 km resolution (left panel), and 60 km resolution (right panel). 



ing with ^ 2000 km s~^ converges to within a factor o f 2-3 (the 
same holds in case of the 120 km model obtained by Guzman 
1200 9). In addition, this result is less sensitive to the width of the 
exclusion cone around the symmetry axis. 

It should be noted that in addition to the variations intro- 
duced by the numerical resolution, and the ambiguity of ana- 
lyzing the simulations near the poles, there are two further un- 
certainties in the models that concern the nickel velocities and 
the initial nickel distribution. First, our axisymmetric simula- 
tions tend to overestima te the d rag coefficient of true three - 
dimensional clumps (K ane et al.l [2000; Ham mer et al.l |200^ . 
And second, the abundance ratios of different iron-group nuclei 
in neutrino-driven supernova models depend very sensitively on 
the neutrino luminosities, and cannot be calculated accurately 
with a grey neutrino transport scheme, as e.g. the one that was 
employed in calculating model b23a (see Paper II). In other 
words, the ^^Ni abundance in model b23a may be lower than in 
reality, favoring nuclei like ^^Fe, ^^Fe, etc., instead. In Paper II 
we have tried to compensate for this uncertainty by considering 
only the total mass fraction of iron group nuclei in the analy- 
sis of velocity distributions. We have proceeded similarly here - 
lumping together different iron group nuclei into what we called 
the "nickel" abundance - since the total abundance of the iron- 
group is expected to be a more reliable quantity than its indi- 
vidual abundances, and can thus serve as an upper limit to the 
^^Ni abundance. A more detailed investigation of these issues is 
required, though, and will be given in future work. 



5. Conclusions 

We presented the results and detailed analysis of a series of ax- 
isymmetric hydrodynamic simulations of the first week of the 
evolution of a non-spherical core-collapse supernova. Different 
from our previous work, our computations were performed in a 
single domain using cylindrical coordinates. We implemented a 
workload-constrained mesh adaption strategy that allowed us to 
complete the simulations given limited computational resources, 
and to avoid the cumbersome periodic mesh remappin g used in 
past studies (iKifonidis et al.ll2006l: [Couch et al.1 12009). We ob- 
tained a series of models at progressively higher mesh resolu- 
tion and provided insight into the numerical convergence of our 
simulations. 

Our simulations are the first to follow a S ASI-dominated ex- 
plosion from shortly after its initiation into the homologous ex- 
pansion and early SNR phase. Especially during the first several 



hours of this evolution, the ejecta are characterized by complex 
interactions between Rayleigh-Taylor, Richtmyer-Meshkov, and 
Kelvin-Helmholtz instabilities, which produce an extensive mix- 
ing and outward penetration of stellar layers enriched in heavy 
elements. The global asphericity of the supernova shock is, 
moreover, essential for triggering a deep inward penetration of 
hydrogen and helium into the central ejecta regions. For the par- 
ticular explosion model and progenitor that we have studied, the 
ejecta become homologous approximately one day after the ex- 
plosion. 

A very important result of our present simulations is that the 
2D SASI instability, which acts during the explosion launching 
phase, not only determines the structure of our model around that 
time, but leaves also a strong, large-scale imprint in the ejecta in 
the form of a significant lateral velocity gradient (Fig.[TT]), which 
affects the evolution for minutes to hours later (i.e. for several 
RT-growth time scales). This means that although the radiative 
driving of the explosion essentially ends around 1 second after 
the core bounce, it is actually incorrect to consider the late-time 
evolution as a phase which is completely detached from the his- 
tory of the explosion. 

This strong lateral expansion of SASI explosion models 
sheds new light onto the formation of the polar outflows which 
have been observed in previous axisymmetric simulations of 
late-time supernova evolution. It had actually been mentioned in 
Paper II that the appearance of these outflows was likely not due 
to a single numerical problem, but rather due to the combination 
of the restriction of the degrees of freedom of the flow, the use of 
reflecting boundaries at the axis of symmetry, and the presence 
of a coordinate singularity, and hence the possible occurrence of 
non-negligible discretization errors at this axis. 

A crucial finding of our present work is that even this quite 
differentiated view for the formation of the polar outflows is 
incomplete. It neglects the fact that also a physical cause is 
involved, namely the strong late-time lateral expansion of 2D 
SASI explosion models. Our simulations indicate that this is ac- 
tually the dominant effect for the formation of such pronounced 
flows near the symmetry axis. The strong lateral expansion away 
from the equatorial plane and toward the poles has, moreover, 
two very important consequences for the observational outcome 
of the models: it results in an inevitable advection of high- 
velocity nickel-rich material from moderate latitudes toward the 
poles, and it contributes to ultimately shape the ejecta into the 
form of a prolate ellipsoid. 
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In the present model, the highest velocity nickel on our grid, 
which moves with speeds close to 4000 km s~^ and is consistent 
with the data of SN 1987 A, is advected within only 300 sec- 
onds after core bounce over half a quadrant of the computational 
domain. In other words, although this nickel-rich material is ini- 
tially located far away from the axis, by the end of the simulation 
it has ended up close to the poles. If confirmed in future three- 
dimensional simulations, this eff'ect of the SASI might actually 
explain the asymmetric nickel lines of SN 19 87 A. 

Yet, since the accuracy of the solution is largely unknown 
near the poles, there is significant ambiguity in what should be 
accounted for in these regions when computing the maximum 
nickel velocities of our 2D simulations. Given that a physical 
eff'ect is involved in transporting material to the poles, it is for 
instance not at all clear whether the use of an exclusion cone 
at the poles - which was customary in earlier simulations - is 
indeed justified. It is also unclear how large such an exclusion 
cone should be made. This ambiguity is impossible to overcome 
in axisymmetric models. 

The medium (i.e. 60 km, and 30 km) resolution models and 
partly also the high (i.e. 15 km) resolution model that we pre- 
sented here, yield results which are in very g ood agree ment with 
the simulations reported in Paper II (.Kifonidis et al.l 12006 ), al- 
though these latter simulations were obtained with a diff'erent 
code, and a diff'erent computational strategy, and were followed 
to only 20 000 seconds after core bounce. In particular we found 
that the amount of mixing of both light and heavy elements, and 
the maximum nickel velocities in both studies agree well. 

The extent of the mixing of heavy elements in our high- 
resolution simulations diff'ers qualitatively from models which 
are less well r esolved (as e.g. the 240 km model obtained by 
lGuzmanll2009l) . At high resolutions, the mass distribution dis- 
plays a pronounced hump with velocities < 1000 km s~^ and 
a long tail extending out to > 3200 km s"^ The hump region 
shows signs of convergence while quantitative diff'erences are 
seen in the tail. At the end of the simulations, the highest nickel 
velocities vary between 3200 km s~^and ~ 4000 km s~\ while 
the mass of the fastest moving nickel varies by at least a factor 2- 
3 between our best resolved models. These results are sensitive 
as to whether the polar regions are included in the analysis or 
not. This shows that due to the strong non-linearity of the prob- 
lem, and the ambiguity in analyzing two-dimensional simula- 
tions near the axis of symmetry, strict numerical convergence is 
difficult to achieve. It also implies that present three-dimensional 
simulations, with their much lower resolution, must be viewed 
as being far from resolved, and that their conclusions must be 
verified by a proper numerical convergence analysis in three di- 
mensions. 

The second important effect of the S ASI-induced lateral ex- 
pansion is that it contributes significantly to the strong prolate 
deformation of the ejecta, as it is observed at the end of our 
simulations (i.e. at ^ = 7 days). This ejecta deformation may 
be considered final, becaus e the expans ion has long become ho- 
mologous by that time. As iK jaer et aP (|2010) have pointed out 
in their recent study, it is moreover in very good agreement with 
the ejecta morphology of SN 19 87 A, making the assumption of 
a "jet-induced" explosion unnecessary. 

Given the importance of the SASI-induced late-time lateral 
expansion for both SNR morphology and the distribution of 
heavy elements, systematic future studies are required to inves- 
tigate how it depends on the dominant SASI mode when the 
early explosion phase ends. In the present simulations it is the 
/ = 2 (quadrupolar) mode of the SASI which ultimately became 
dominant, but other modes are clearly possible in 2D simula- 



tions (IScheck et al.l l2006h . while little is known in the three- 
dimensional case. Since a thorough survey of SASI mode out- 
comes cannot be completed in two dimensions and considering 
furthermore the problems that the analysis of 2D simulations 
poses with regard to the determination of the maximum nickel 
velocities, 3D high-resolution simulations should be performed 
as soon as corresponding computational resources will allow for 
that. Such three-dimensional simulations will need to employ 
singularity-free grids that cover the entire sphere, in order to be 
free of the aforementioned problems, and hence of optimal use 
in both early and long-time supernova modeling. 
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